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SECTION  I:  GENERAL  PLASMA  THEORY  AND  SIMULATION 

A.  Linear  Magnecized  Plasma  Response  to  an 
Oblique  Electrostatic  Wave 

William  S.  Lawson 

Abstract 

The  linear  response  of  a  spatially  periodic  magnetized  Vlasov 
plasma  distribution  function  is  computed  to  second  order  in  the  elec¬ 
tric  field.  The  results  for  a  specific  electric  field  are  then  compared 
with  the  results  of  computer  simulation  for  different  amplitudes  of  the 
electric  field.  Both  trapping  and  resonant  heating  are  observed,  and 
both  appear  to  contribute  (for  the  chosen  parameters)  to  limiting  the 
validity  of  linear  theory  at  larger  electric  field  amplitudes. 

This  is  to  be  issued  as  ERL  Memo  No.  UCB/ERL  M86/94,  Dec.  15,  1986. 

The  abstract  follows. 

B.  Limits  of  Linear  Response  of  a  Vlasov  Distribution 

William  S.  Lawson 

The  abstract  for  this  appeared  in  the  previous  QPR.  The  report  is 
being  issued  as  ERL  Memo  No.  UCB/ERL  M86/44,  June  3,  1986. 

C.  Saturation  and  Post-Saturation  Behavior  of  the 
Alfven  Ion-Cyclotron  Instability:  A  Simulation  Study 

Niels  F.  Otani 

This  is  to  appear  as  ERL  Memorandum  No.  M86/15,  Feb.  19,  1986.  The 
abstract  follows. 

D.  Stabilizing  Effects  of  Finite-Amplitude 
RF  Waves  on  the  Interchange  Instability 

Niels  F.  Otani 

This  is  to  appear  as  ERL  Memorandum  No.  UCB/ERL  M86/18,  Feb.  28,  1986. 
The  abstract  follows. 

E.  Effect  of  Large-Amplitude  Perpendicularly-Propagating 
RF  Waves  on  the  Interchange  Instability 

Neils  F.  Otani  and  Bruce  I.  Cohen 

This  is  to  appear  as  ERL  Memorandum  No.  UCB/ERL  M86/48,  July  8  1986 

The  abstract  follows. 


F.  Numerical  Simulations  of  Turbulent  Trapping  in  the  Weak 
Beam-Plasma  Instability:  Study  of  the  Growth  Rates 

Dr.  Kim  Theilhaber 

This  work  will  be  issued  in  an  ERL  Memorandum,  "Numerical  Simulations 
of  Turbulent  Trapping  in  the  Weak  Beam-Plasma  Instability,"  by  K.  Theilhaber 
ERL  Memo  M86/50,  June  5,  1986.  The  abstract  follows. 


Saturation  and  Poit-Satnration  Behavior  of  the 
Alfven  Ion-Cyclotron  Instability:  A  Simulation  Study 


Niels  F.  Otani 

Electronics  Research  Laboratory 
University  of  California 
Berkeley,  California  94720 


Two  computer  codes  based  on  different  algorithms  are  employed  in  the  study 
of  the  Alfven  ion-cyclotron  instability.  Features  of  quasilinear  diffusion  are  clearly 
exhibited  and  identically  modeled  by  the  two  codes.  Ion  trapping  is  also  observed 
and  appears  to  play  an  important  role  in  the  saturation  process.  Relaxation 
of  the  anisotropy  in  the  ion  temperature  and  transverse  wave  energy  transfer 
from  intermediate  to  long  wavelength  modes  are  observed  to  proceed  more  slowly 
following  instability  saturation  when  the  possibility  of  nonlinear  E  x  B  electrons 
is  allowed,  and  more  slowly  still  when  longitudinal  ion  sound  modes  are  present. 
Fluctuations  in  the  short  wavelength  modes  are  shown  to  be  a  discrete  particle 
effect,  with  levels  agreeing  with  predictions  calculated  from  a  test-particle  model. 


Stabilizing  Effects  of  Finite-Amplitude  RF  Waves 
on  the  Interchange  Instability 


Niels  F.  Otani 

Electronics  Research  Laboratory 
University  of  California 
Berkeley,  California  94720 


The  effects  of  the  presence  of  a  finite-amplitude  RF  wave  on  the  interchange 
instability  are  studied  theoretically  in  a  three-dimensional  quasineutral  Darwin 
two-fluid  system.  An  exact  energy  conservation  law  is  developed  for  the  system 
equations.  The  perturbed  dynamical  quantities  are  expanded  to  second  order 
in  the  fluid  displacements  and  it  is  demonstrated  that  a  Lagrangian  exists  for 
the  system.  The  interchange  perturbation  energy  is  calculated  and  terms  associ¬ 
ated  with  the  perturbed  magnetic  potential  energy  are  analyzed.  For  the  terms 
studied,  it  is  found  that  no  stabilization  of  the  incompressible  interchange  mode 
occurs  unless  the  RF  wave  has  a  parallel-propagating  component,  and  that  cou¬ 
pling  to  perpendicular-propagating  sidebands  does  not  influence  the  stability  of 
perpendicularly-propagating  interchange  modes.  The  stabilizing  influence  of  these 
terms  also  does  not  depend  directly  on  the  presence  of  RF  field  gradients.  The 
analysis  suggests  a  physical  mechanism  by  which  RF  stabilization  can  occur — the 
RF  wave  stabilizes  the  interchange  mode  by  forcing  the  mode  to  bend  magnetic 
field  lines. 


Effect  of  Large-Amplitude  Perpendicularly-Propagating 
RF  Waves  on  the  Interchange  Instability 

Niela  F.  Otani* 

Electronics  Research  Laboratory 
University  of  California 
Berkeley,  California  94720 

Bruce  I.  Cohen 

Lawrence  Livermore  National  Laboratory 
University  of  California 
Livermore,  California  94550 


Results  are  presented  from  hybrid  2-d  quasineutral  Darwin  simulations  of  the 
interchange  instability  in  the  presence  of  a  large  RF  wave  in  the  ion-cyclotron  fre¬ 
quency  range.  The  simulation  models  the  plane  perpendicular  to  the  background 
magnetic  field  using  cold,  particle  ions  and  a  cold  E  x  B  electron  Quid.  Related 
theory  is  also  discussed.  Fluid  equations  appropriate  to  our  simulation  model  are 
derived  and  their  properties  demonstrated  and  compared  to  simulation.  A  method 
for  solving  for  the  RF-modified  growth  rates  from  the  fluid  equations  is  described. 
It  is  generally  expected  that  the  current  component  associated  with  the  mean,  RF- 
induced  ion  drift  is  capable  of  influencing  the  stability  of  the  interchange  mode; 
however,  no  modification  of  the  mean  ion  drift  is  observed  in  simulations  in  which 
RF  is  present.  Instead,  in  both  the  theory  and  simulation,  an  electron  RF-field 
oscillation  current  dominates  the  modification  to  the  gravitational  current.  As 
a  result,  even  in  the  presence  of  large  RF  fields,  ( Bur /Bo  =  15%)  only  modest 
corrections  to  the  interchange  growth  rates  are  observed.  The  effect  is  stabilizing 
for  kLn  <  0. 8-0.9,  apparently  for  both  signs  of  the  square-electric-field  gradient, 
and  is  destabilizing  for  larger  values  of  kLn,  although  the  credibility  of  the  sim¬ 
ulation  begins  to  become  suspect  here.  Fractional  reduction  of  the  interchange 
growth  rate  is  observed  to  be  quadratically  dependent  on  the  RF  wave  amplitude, 
independent  of  ion-cyclotron  resonant  effects,  and  proportional  to  V B\f/V Bq, 
consistent  with  an  eikonal  theory  developed  for  the  study  of  stabilizing  effects  on 
perpendicularly-propagating  fast  Alfven  waves.  The  results  also  suggest  that  ad¬ 
ditional  gradient-independent  stabilizing  effects  may  be  operative  when  kLn  —  1. 
Finally,  it  is  also  observed  that,  while  the  RF  wave  has  little  effect  on  the  inter¬ 
change  instability,  the  interchange  mode  strongly  affects  the  RF  wave,  damping 
it  significantly  as  the  mode  saturates. 


*  Present  address:  Courant  Institute  of  Mathematical  Sciences,  New  York  University, 
New  York,  New  York  10012. 


NUMERICAL  SIMULATIONS  OF  TURBULENT  TRAPPING  IN  THE 
WEAK  BEAM-PLASMA  INSTABILITY:  STUDY  OF  THE  GROWTH  RATES 


K.  Theilhaber 

Plasma  Theory  and  Simulation  Group 
Electronics  Research  Lab 
University  of  California 
Berkeley,  CA  94720  USA 


G.  Laval  and  D.  Pesme 
Centre  de  Physique  Theorique 
Ecole  Polytechnique 
91128  Palaiseau,  France 


ABSTRACT 


Numerical  simulations  of  the  weak  beam-plasma  instability 
were  done  in  the  turbulent  regime  where  small-scale  trapping  is 
a  dominant  feature  of  the  instability,  a  regime  with  behavior 
not  predicted  by ^qu^s^linear  theory.  It  occurs  when  the  trapping 
frequency  v  *  (k  D)  is  larger  than  the  growth  rate  of  the 
instability.  The  results  of  the  simulations  were  compared  with 
those  of  a  specific  model  of  the  turbulence,  the  so-called 
"turbulent  trapping"  model,  which  predicts  a  growth  rate  well 
enhanced  over  the  quasilinear  value.  It  was  found  that  while 
growth  rates  were  enhanced  over  the  quasilinear  values,  the 
enhancements  observed  are  smaller  than  expected  from  the  quanti¬ 
tative  predictions  of  the  model.  Further  work  is  necessary  to 
determine  whether  this  discrepancy  is  a  failing  of  the  turbulent 
trapping  model,  or  the  result  of  the  numerical  limitations  of  our 
computational  scheme. 
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SECTION  II:  PLASMA-WALL  PHYSICS,  THEORY  AND  SIMULATION 


A.  Planar  Magnetron  Discharge  Simulations 
Perry  Gray  (Prof.  C.  K.  Birdsall) 

A  one-dimensional  particle  plasma  simulation  code,  following  the 
self-consistent  motion  of  thousands  of  particles,  has  been  developed  to 
complement  theory  and  experiment  in  studies  of  plasma-sheath  and  sheath- 
wall  interactions,  for  a  variety  of  applications  with  respect  to  plasma- 
assisted  materials  processing  problems. 

We  have  implemented  volume  generation  of  ions  and  electrons  (in 
contrast  to  the  surface  emission  model  described  earlier) .  The  ionization 
rate  can  be  made  a  function  of  some  local  variable  (e.g.,  local  electric 
field  or  electron  density).  We  have  also  added  transport  diagnostics  in 
the  form  of  fluxes  (particle,  kinetic  energy,  and  heat)  and  temperature 
as  functions  of  x  and  secondary  emission  from  the  cathode  surface.  We  will 
add  to  the  model  the  effects  of  charge-neutral  collisions.  The  bounce  mo¬ 
tion  of  the  particles,  from  mirror  to  mirror,  may  be  added.  With  these 
additions,  we  believe  that  much  of  the  essential  physics  of  the  crossed 
field  discharges  will  be  modeled.  Initial  simulation  studies  of  the  dis¬ 
charge  have  been  performed. 

The  first  simulation  studies  performed  use  the  one-dimensional  (1-d, 
3-v)  model  of  the  planar  magnetron  shown  in  Figure  1.  The  ignorable  dir¬ 
ections  are  y  and  z.  The  dimension  of  the  simulation  Is  a  cut  along  x, 
the  axis  of  symmetry,  so  that  the  magnetic  field  is  perpendicular  to  the 
direction  of  simulation.  We  keep  track  of  velocities  in  all  three  direc¬ 
tions  so  that  gyromotion  of  the  ions  and  electrons  is  included.  The 
magnetic  field  B  (x)  is  assumed  to  be  due  to  a  linear  dipole,  offset  with 
respect  to  the  right  boundary  of  the  simulation  as  shown,  or  to  be  uniform. 
The  left  boundary  of  the  simulation  is  at  ground  (<J  =  0)  ,  while  the  right 
boundary  has  a  fixed  d.c.  bias  (and  may  also  be  a.c.  driven). 

Typical  laboratory  parameters  are: 

kT  «  5eV,kT.  *  0.025eV  (T.<c  T  ,  essentially  cold  ions) 

e  1  l  e 

~  ,n12/  3 

n.  =n  =  10  /cm 

l  e 

B  =  lkG 
z 

which  produces: 

( ju  ~  6x10^  sec  * ,  CO  -sr?  1.8x10^  sec  * 

pe  ce 

8  \ 

v  ^  10  cm/sec,  A.,  1.6x10  cm 

te  De 


0. 


Typical  initial  results  for  B  uniform  in  (x)  are 
shown  in  Figure  2.  The  left-hand  electrode  is 
grounded  (Sf  m  0)  and  the  right-hand  electrode  is  held 
at  e0/ lcTe  ■  -20.  (Typical  laboratory  discharges  have 
kTe  ^  5eV,  so  this  is  a  bias  of  -100  volts.)  The  ions 
are  considered  cold.  (Typically  kT^  as  0.025eV  so 
Ti^<  .)  Ue  use  a  mass  ratio  of  m^/me  ■  40  so  as  to 
be  able  to  follow  both  electrons  and  ions  conveniently. 

The  results  shown  are  at  time  t  ft;  31  (after  4000 
time  steps)  with  roughly  2000  each  of  electrons  and 
ions.  The  potential  0(x)  remains  near  0  over  most 
at  the  space,  then  starting  at  x  8  dips  to 
e  0(L)/kT_  ■  -20.  The  electrons  are  repelled  by 
the  positive  electric  field  in  this  space  x>8,  as 
observed  in  the  vx  -  x  phase  space  plot,  never  reach¬ 
ing  beyond  x«  8.  On  the  other  hand,  the  ions,  be¬ 
ginning  at  x  ■  8,  are  accelerated  by  the  E  ,  hit¬ 
ting  the  right  hand  plate  with  high  velocity,  well 
beyond  vfj  (accelerated  by  close  to  the  100  volts) . 
One  quantitative  feature  is  that  the  curvature  of 
0(x)  is  negative  so  that  the  power  law  is 
60(x)/v  (x-8)n  with  n>l;  similarly  the  curvature 
of  the  ion  acceleration  is  like  v^  (x)-v(x-8)m  with 
m  <  1.  The  temptation  is  to  say  that  n  -  4/3,  a 
Child's  Law  region,  making  m  ■  2/3  as  the  ions  are 
accelerated  nearly  from  rest.  That  is  to  say,  the 
ions  free-fall  near  the  right  hand  electrode,  in¬ 
fluenced  by  ,  as  if  Bz  were  zero.  The  correla¬ 
tion  is  good  with  experiment  in  that  a  dark  space  is 
observed  next  to  the  negative  plate  (where  the  elec¬ 
tron  population  is  quite  low)  and  its  length  is 
associated  with  the  sheath  thickness,  and  not  the 
ion  gyroradius .  The  electrons  have  gyroradlus  of 
ag-v  /uice*l/8  and  the  thermal  ions  have  a  gyroradius 
of  aj  -  vti/ujcl-  0.005/(8/40)  -  0.025,  both  fairly 
well  magnetized.  Hence,  how  do  the  ions  free  fall, 
if  they  are  so  tightly  bound  to  Bz  (gyroradius  small 
compared  with  L)  ?  Note  that  at  x%8,  the  cold  ions 
have  become  accelerated  a  little,  to  about  ^  0 . 1 
which  gives  a  much  larger  a*  (by  the  factor 
0.1/0.005  ■  20)  of  about  0.5;  near  x  -  9  the  ions  are 
at  v^  a>  0.5  with  again  much  larger  a^  ss  2.5  and  x  -  10, 
ions  are  near  v^rs  1.0  with  a^  -  5  (or,  gyrodiameter 

■  L,  diode  spacing).  Hence,  aj^  increases  dramatically , 
so  that  free  fall  takes  place.  The  transport  of  ions 
needs  no  collisions  to  take  place. 
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.  Planar  Magnetron  Discharges:  Mechanisms  for  Instability,  and 
Numerical  Study  of  Equilibrium  -  K.  Theilh&ber. 


1  Introduction 

In  the  following  report,  we  shall  present  a  two-dimensional  simulation  model 
for  the  planar  magnetron  which  embodies  a  plausible  mechanism  for  insta¬ 
bility  and  fluctuating  behavior  of  the  magnetron  discharge.  In  this  model, 
the  magnetic  field  is  perpendicular  to  the  plane  of  motion  of  the  ions  and 
the  electrons,  in  the  so-called  "crossed- field”  configuration.  This  article  is 
meant  as  a  preliminary  study  which  should  help  put  the  problem  in  an  ap¬ 
propriate  perspective  for  modelling;  it  includes  a  short  review  of  previous 
work  on  instabilities  occuring  in  cross-field  configurations,  some  analytic 
work  specific  to  our  planar  magnetron  model,  and  numerical  work  on  the 
equilibrium  of  this  model.  Our  long  term  goal  will  be  to  simulate  the  entire 
nonlinear  evolution  of  the  two-dimensional  instabilities  which  can  occur  in 
the  planar  magnetron  configuration. 

2  Review  of  Existing  Crossed-Field  Theory 

In  the  past,  crossed-field  instabilities  have  been  examined  in  at  least  two 
configurations,  whose  similarities  with  the  planar  magnetron  geometry  lead 
us  to  believe  that  analogous  behavior  will  occur  in  the  magnetron  discharge 
as  well.  In  both  cases,  the  source  of  free  energy  resides  in  the  shear  of  the  ve¬ 
locity  field  of  the  particles,  a  shear  which  is  induced  by  nonuniform  E  x  B 
drifting  perpendicular  to  the  magnetic  field.  The  two  configurations  we 
shall  mention  embody  extreme  examples  of  crossed-field  instability,  and  we 
expect  that  the  magnetron  model  will  exhibit  an  instability  of  an  interme¬ 
diate  type,  which  partakes  of  both  extremes. 


2.1  The  Diocotron  Instability 

This  is  an  instability  which  occurs  in  a  single-species,  non-uniform,  non¬ 
neutral  plasma,  and  which  has  been  studied  in  connection  with  the  behavior 
of  electronic  devices.  The  model  of  this  instability  allows  for  small  space 
scales  (thin  sheets  of  charge,  etc),  and  the  presence  of  material  boundaries. 

The  discovery  of  the  Diocotron  instability  goes  back  to  work  by  O. 
Buneman[l].  More  recent  work  includes  that  of  Knauer[2]  and  Byers[3]. 
The  simplest  configuration,  which  we  shall  discuss  here,  is  that  of  a  very 
thin  charged  sheet  of  electrons,  drifting  perpendicular  to  a  magnetic  field, 
a  situation  treated  analytically  by  R.  Gould  in  1957  [4],  who  also  provided 
a  physical  explanation  for  the  instability.  The  dispersion  relation  of  such  a 
system  is  found  to  be: 


U>  =  fcyUO  ±  I 


kyO0 

2-BqCo 


(1) 


whre  u0  is  the  mean  velocity  of  the  sheet,  o0  its  charge  per  unit  area,  B0 
the  magnetic  field  strength,  and  ky  the  wave-number  of  the  perturbation 
parallel  to  the  sheet  and  perpendicular  to  B.  If  w  is  the  sheet  thickness, 
and  no  the  mean  electron  density  in  the  sheet,  then  the  growth-rate  can  be 
written: 


1  =  \  l*.»l  (2) 

l  Wgf 

where  u*e  =  n0e2/m,f0l  u>ei  =  eB0/m„  and  with  the  thin-sheet  condition, 
|kyti>|  1.  For  |ky«>|  ~  1,  Eq.(2)  can  be  taken  as  qualitatively  correct.  If 
we  assume  that  wp«  ~  uet,  a  situation  possible  for  the  magnetron  discharge, 
Eq.(2),  which  models  a  very  extreme  configuration,  predicts  a  large  growth 
rate,  ^  ~ 


2.2  The  MHD  Kelvin-Helmholtz  Instability 

This  is  an  instability  which  occurs  in  a  two-species,  almost  neutral  plasma, 
and  which  has  been  studied  in  connection  with  space-science  problems. 
This  model  calls  for  large  scale  lengths  and  boundary  conditions  at  infinity 
(or  very  far  away  from  the  shear  layer). 


In  form,  the  MHD  Kelvin-Helmholtz  instability  is  closer  to  a  pure  hy¬ 
drodynamic  instability  than  the  Diocotron  instability.  It  has  been  studied 
by  Chandrasekhar [5],  who  has  given  explicit  analytic  expressions  for  the 
growth  rate,  under  the  assumption  of  simple  velocity  shear  profiles.  More 
recent  work  can  be  found  in  the  articles  by  Miura  and  Pritchett  [6]  and 
Pritchett  and  Coroniti[7j.  The  first  reference  provides  a  fluid  treatment  of 
the  Kelvin-Helmholtz  instability,  with  the  growth  rates  found  through  a 
numerical  integration  of  the  fluid  equations.  The  second  reference  extends 
the  fluid  analysis  by  presenting  the  results  of  a  particle  simulation,  against 
which  the  predictions  of  the  fluid  picture  can  be  checked. 

In  both  of  [6,7],  the  conditions  assumed  for  the  MHD  model  are  (Fig.(l)): 

1.  an  almost  neutral  plasma  (in  sharp  contrast  to  the  Diocotron  insta¬ 
bility). 

2.  infinitely  magnetized  electrons. 

3.  the  dynamics  are  determined  by  the  ions  (as  opposed  to  the  Diocotron 
instability,  where  electrons,  the  only  species  present,  determine  the 
growth  rates). 

4.  the  ordering  of  the  scale  lengths  is: 

ptC  pi  <  o  <  A„,  (3) 

where  (see  Fig.(l)),  pt  and  p{  are  the  Larmor  radii,  a  the  width  of  the 
shear  layer,  and  A„  the  wavelength  of  the  perturbation  parallel  to  the 
velocity  direction. 

With  the  above  conditions,  and  with  the  velocity  shear  determined  by: 

vo„(x)  =  - E0t{x)/B0 ,  (4) 

where  E0t  is  the  equilibrium  electric  field,  the  analysis  of  the  fluid  equations 
leads  to  a  dispersion  relation  in  which  the  frequency  lj  is  obtained  as  a 
solution  of  the  eigenvalue  problem: 


where  6viz  is  the  velocity  perturbation,  with  6t/,-,(±oo)  =  0.  Corresponding 
to  each  value  of  a/,  Eq.(5)  predicts  an  eigenfunction  6t/,*(x). 

Miura  and  Pritchett[6]  obtained  the  growth-rate  for  the  shear  profile: 

Voy(x)  =  v0  tanh (x/a),  (6) 

where  a  determines  the  width  of  the  shear  layer  and  Vo  its  magnitude.  In 
fact,  the  authors  of  [6]  did  not  resort  to  a  direct  solution  of  Eq.(5),  but 
rather  to  a  numerical  integration  in  time  of  the  system  of  linearized  fluid 
equations:  the  resulting  profile  for  the  growth  rate  is  shown  in  Fig.(2). 
The  growth-rate  has  a  maximum  for  kva  ss  0.5,  for  which  Av  as  4*a  as 
12a.  As  the  eigenfunctions  extend  over  a  distance  in  x  comparable  to  Ay, 
the  maximum  growth  rate  occurs  for  velocity  perturbations  extending  over 
many  times  the  shear  length.  In  its  nonlinear  stages,  the  instability  develops 
vortices  which  are  also  considerably  larger  than  the  dimensions  of  the  shear 
layer. 

A  surprising  result  of  the  particle  simulations  of  Pritchett  and  Coroniti[7] 
is  that  they  show  that  the  fluid  calculations  outlined  above  are  valid  even 
when  the  ion  gyroradius  is  comparable  in  size  to  the  width  of  the  shear 
layer,  a  ~  pi  (though  not  when  it  is  smaller,  a  C  p*).  This  is  a  situation 
where  not  only  one  expects  finite-Larmor-radius  effects  to  be  large,  but  in 
which  it  is  extremely  hard  to  model  these  effects,  short  of  writing  the  full 
kinetic  equations.  Thus,  the  particle  simulations  encourage  us  to  think  that 
the  fluid  formulation  is  a  more  powerful  tool  than  a  conservative  estimate 
would  lead  us  to  believe,  and  that  it  can  be  used  for  a  first  estimate  of 
growth  rates  in  similar  configurations. 

3  The  Magnetron  Model 

We  shall  now  outline  a  "first”  magnetron  model.  Our  approach  in  defining 
the  model  is  to  "streamline"  it  as  much  as  possible,  so  that  in  its  final  form 


it  will  exhibit  the  crossed-field  instability  in  the  simplest  form  attainable, 
while  retaining  the  basic  physics  of  the  process. 

The  first  problem  in  designing  the  model  is  to  determine  the  equilibrium 
conditions.  The  creation  and  sustenance  of  the  magnetron  discharge  are 
thought  to  strongly  relie  on  collisional  processes,  such  as  neutral  ionization 
and  secondary  electron  emission,  which  depend  not  only  on  plasma  prop¬ 
erties,  but  on  the  properties  of  the  ambient  neutral  gas,  and  of  the  plate 
surfaces.  The  discharge  mechanism  is  illustrated,  in  a  schematic  and  tenta¬ 
tive  form,  in  Fig.(3).  We  believe  that  trying  to  include  all  of  these  primarily 
collisional  effects  in  the  model,  is  premature  at  our  stage  of  understanding, 
and  we  thus  choose  to  isolate  the  model  from  the  detailed  nature  of  the 
source  mechanisms,  by  postulating  instead  ad  hoe  sources  (planar  or  dis¬ 
tributed)  for  the  ions  and  electrons.  The  simplified  geometry  is  shown  in 
Fig. (3).  As  we  shall  see,  even  the  choice  of  these  given  sources  is  not  sim¬ 
ple,  as  the  instability  mechanism  is  very  dependent  on  the  precise  nature  of 
the  sources  (injection  rates,  distribution  in  space).  We  are  thus  confronted 
with  the  problem  of  choosing  the  most  plausible  source  condition,  in  the 
absence  of  a  detailed  understanding  of  the  magnetron  discharge. 

A  further  simplification  of  the  model  is  to  concentrate  on  the  sheath 
region  represented  in  Fig. (3),  across  which  occurs  almost  all  of  the  volt¬ 
age  drop.  Its  representation  in  the  simulation  model  is  shown  in  Fig. (4). 
Further  assumptions  of  the  simulation  model  are: 

1.  we  model  only  the  cross-field  dynamics  of  the  sheath  (in  the  plane 
perpendicular  to  the  magnetic  field). 

2.  the  sheath  plasma  is  a  two-species,  non-neutral,  highly  nonuniform 
plasma. 

3.  in  what  follows,  the  ion  source  is  planar:  the  ions  are  emitted  from 
the  right-hand,  high-potential  plate  as  if  from  a  "plasma  reservoir”. 
The  ions  are  emitted  almost  cold,  with  T,  <c  Te. 

4.  temporarily,  we  shall  treat  the  electrons  as  initially  distributed  in  the 
simulation  region,  with  no  subsequent  electron  injection. 

5.  the  electron  gyroradius  is  much  smaller  than  the  sheath  thickness  L: 
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6.  the  ions’  radius  of  curvature  in  falling  to  the  collector  plate  is  larger 
or  of  the  same  order  as  the  sheath  thickness  L. 

The  model  outlined  above  shares  features  of  the  models  for  the  Dio- 
cotron  and  the  MHD  Kelvin-Helmholtz  instabilities:  it  allows  for  a  highly 
non-neutral,  bounded  plasma,  but  one  in  which  the  ion  dynamics  play  an 
essential  role. 

Our  approach  will  be  to  analyse  the  equilibrium  properties  of  the  sheath, 
by  first  removing  all  electrons  from  the  system.  This  leads  to  an  equilibrium 
which  is  nothing  but  a  Child-Langmuir  sheath,  and  for  which  the  numerical 
results  can  easily  be  checked  against  the  well-known  analytic  results.  We 
shall  then  consider  analytically  the  effect  of  a  nonuniform  distribution  of 
electrons. 


4  Equilibrium  without  Electrons:  the  Ion 
Child-Langmuir  Sheath 

We  shall  outline  the  analysis  for  the  Child-Langmuir  sheath,  in  the  ab¬ 
sence  of  a  magnetic  field,  and  assuming  almost  cold  ions.  We  assume  that 
the  plasma  "reservoir”  to  the  right  of  the  high-potential  plate  injects  ions 
according  to  the  flux: 


r,(L)  =  —0  ni0vti ,  0  =  1/(2tt)1/j  «  0.4,  (7) 

With  the  flux  determined  by  Eq.(7),  the  local  ion  density  is  given  by: 

0niOvtt 


n<(x)  = 


(8) 


(2  «/m,)1/J(V  —  <£(x))*/2’ 
where  V  is  the  potential  at  the  emiiting  plate.  We  add  Poisson’s  equation: 

,  _  _  e_ 

dx'*  ton" 


and  define: 


(9) 


£  =  z/Zr,  \l>  =  <t>/V, 


(10) 


vo  =  (2e/mj)1/sVl/a,  (11) 

Then,  for  a  monotone  potential  profile,  we  obtain  the  solution  0(£)  as  the 
solution  of  the  implicit  equation: 

(*+ r?)17’  (*  -  rO  -  (**•<« + W'  -  j«?)  =  I*’  «■ 

(12) 

where  F  =  (1  —  0)1//a,  and  where  the  parameter  K  determines  the  depen¬ 
dence  of  the  solution  on  the  ion  parameters: 

Ks 

Co  V  Vo 

The  integration  constant  t\  is  determined  by  the  condition: 

(^  +  Ie;)'/’(K_I,;)  +  ie;=  ?*•>,  (i4) 

Now  for  fixed  voltage  and  sheath  length,  the  parameter  K  will  increase  as 
the  injection  rate  increases.  There  is  a  critical  value  of  K,  Kc  =  4/9  « 
0.444,  for  which  Ez  =  0  at  x  =  L,  and  for  which  we  get  the  "traditional” 
Child-Langmuir  sheath  potential  and  electric  field  profiles: 

*(*)  =  V  (l  -  (1  -  x/L)‘>3)  ,  (15) 

E,{z)  =  -~{l-xlL)'l\  (16) 

For  K  >  Kc  (and  vti  v0),  Eqs.(  15,16)  remain  valid:  injection  of  ions  at 
x  —  L  has  "saturated”,  and  the  excess  ions  are  turned  back  to  the  injection 
plate. 
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5  Numerical  Results  for  Equilibrium 

We  performed  numerical  runs  of  the  two-dimensional,  electrostatic  code 
ES2(8],  with  only  ions  present,  with  the  intent  of  checking  numerical  results 
against  Eqs.(l5,16).  In  Figs. (5, 6)  we  show  the  results  for  the  simulation 
with  the  geometry  shown  in  Fig. (4).  This  configuration  is  unmagnetized, 
with  B  =  0.  The  plasma  parameters  are:  t >u  =  0.079,  v0  =  (2 eV/m,)1/2  = 
0.98,  e  =  0.01,  e/m,  =  1/40,  and  n,0  =  303.  For  this  density  of  the 
ion  reservoir,  K  =  Kc  =  4/9.  The  external  circuit  has  V  =  19.5  and 
R  =  L  =  0.  The  agreement  for  the  potential  and  electric  field  profiles  is 
quite  close,  except  near  the  emitting  plate  at  x  =  10,  where  the  assumption 
of  completely  cold  ions  (7<  =  0)  is  not  strictly  valid.  In  particular,  the 
simulation  value  of  the  electric  field  at  the  emiiting  plate  is  nonzero,  with 
EX(10)  «  —0.55.  In  fact,  in  Fig.(7),  we  show  that  as  we  reduce  the  injected 
ion  temperature,  the  current  saturates  to  a  value  closer  to  the  ideal,  cold-ion 
saturation. 

We  have  also  studied  the  effect  of  a  nonzero  magnetic  field.  From  simple 
kinematic  arguments,  one  can  argue  for  the  existence  of  a  critical  magnetic 
field  strength  Be,  above  which  all  ion  current  to  the  collector  should  be 
suppressed.  This  field  is  determined  by  the  requirement  that  an  ion  leaving 
the  emitter  plate  with  zero  velocity  will  just  barely  graze  the  collector  plate 
before  being  turned  back  by  the  magnetic  field.  It  is  found  to  be: 


to,  (2 eV /fm)1'2 

Be  e  L  ’  (17) 

This  effect  is  illustrated  in  Fig.  (8),  where  we  show  the  current  as  a  function 
of  increasing  magnetic  field,  with  otherwise  the  same  simulation  parame¬ 
ters  as  in  Fig. (5).  Scatter  plots  of  the  particles  in  phase  space  are  shown  in 
Figs. (9),  showing  how  the  ion  trajectories  turn  over  with  increasing  mag¬ 
netic  field. 
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Putting  In  Electrons  -  Preliminary  The¬ 
ory 


The  electric  fields  predicted  by  Eqs.(15)  or  (16)  are  highly  non-uniform 
across  the  sheath  region.  If  we  think  of  adding  an  electron  population  to 
the  sheath,  perhaps  initially  in  small  enough  numbers  so  as  not  to  notably 
modify  the  ion-sheath  equilibrium,  they  will  acquire  an  Ex  B  motion  which 
will  be  sheared  in  x,  and  this  new  equilibrium  distribution  presents  the  pos¬ 
sibility  of  Kelvin-Helmholtz  instabilities.  In  introducing  electrons  into  the 
model,  we  are  however  faced  with  the  difficulty  of  modeling  their  injection 
into  the  sheath  region.  If  we  attempt  to  inject  them  at  the  high-voltage 
plate,  in  the  same  manner  as  the  ions,  they  will  not  penetrate  further  than 
a  gyroradius  into  the  simulation  region,  leading  to  an  unphysical  electron 
sheet  at  x  =  L.  On  the  other  hand,  if  we  introduce  the  electrons  as  an  ini¬ 
tially  localized  population,  say  in  the  form  of  a  slab,  we  can  expect  a  large 
growth  rate  of  the  order  of  the  "pure”  Diocotron  instability.  Thus  the 
electrons  will  be  rapidly  dispersed,  and  we  suspect  that  dynamic  behavior 
will  cease  as  soon  as  the  elctron  distribution  has  been  smeared  out.  Thus 
it  would  seem  that  to  model  a  fluctuating  steady-state,  with  a  distributed 
electron  population,  will  require  a  distributed  and  continuous  creation  of 
electrons  in  the  sheath  region.  In  what  follows,  we  shall  merely  present 
a  tentative  theory  for  the  initial-value  problem,  in  which  a  small  electron 
population  is  initially  added  to  the  ion-sheath. 

We  assume  that  the  initial  electron  density  distribution  is  of  the  form, 

n0(x)  =  N0  +  (18) 

and  that  the  equilibrium  electric  field  £o(x)  is  known  as  well,  so  that  the 
equilibrium  drift  velocities  of  the  electrons  are  given  by: 

u0v(x)  =  E0{x)/B0,  (19) 

We  shall  use  the  continuity  equation: 

dn 

—  +  V  •  (nv)  =  0,  (20) 
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the  momentum  equation,  which  reduces  to: 


Ext 


Bo  ' 


and  the  Poisson  equation, 


vv  =  ~, 


where  n  refers  to  the  electron  dynamics.  We  assume  for  the  moment  ion 
independent  dynamics. 

If  we  linearize  according  to: 


n  =  n0(i)  +  ni, 


v  =  v0  +  Vi, 


then,  using  the  fact  that  V  •  Vi  =  0,  we  find: 


(dt  +  voydjni  =  -vu-^7> 


vix  —  B\y  /  Bq, 


Ely  —  —dv<f>l, 


d\<t>  i  =  —  nu 


This  system  of  equations  reduces  to: 


where  the  prime  denotes  differentiation  with  respect  to  x.  With  : 
c3Ni/Me(o,  and  with: 


L  uct 
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fl(l)  =  U>  -  kyV  0y(x), 

(31) 

Eq.(29)  takes  the  form: 

(32) 

where  M  is  constant.  Eq.(29)  has  no  obvious  analytical  solutions,  and  it 
would  seem  that  an  approach  analogous  to  that  of  [6],  with  the  numerical 
integration  in  time  of  Eqs. (25-28)  to  determine  the  growth  rate,  might  be 
preferable  to  a  direct  attack  on  Eq.(32).  However,  we  shall  close  this  Section 
by  suggesting  a  variational  approach  to  the  solution  of  Eq.(32).  If  we  write: 


il>(x)  =  <t>{  x)/n(x), 

then  Eq.(32)  takes  the  self-adjoint  form: 

s(n’^)+n(n"-*«’n  +  M)  =  0- 

Multiplying  Eq.(34)  by  ip'  one  then  obtains  a  quadratic  equation: 

Au>*  +  Buj  +  C  =  0, 


where: 
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dx 


Loo  dx  *vvov(i)  (M  -  k„v J„(x)  +  2ksvv0v(x))  , 


-  20 


(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 


Eq.(35)  is  in  fact  a  variational  principle  for  the  eigenfrequency  u>.  One  can 
then  show  that  the  condition  C  <  0  for  all  possible  trial  functions  i>  is  a 
sufficient  condition  for  stability.  However,  the  expression  for  C  is  certainly 
not  negative  definite,  and  we  can  presume  that  Eq.(35)  has  unstable  roots. 

7  Conclusion  and  Future  Work 

We  have  outlined  above  the  framework  for  numerical  and  analytic  (or  semi- 
analytic)  work  on  cross-field  instabilities  in  a  configuration  which  models 
part  of  the  planar  magnetron.  Our  basic  magnetron  model  is  not  yet  com¬ 
plete,  the  remaining  problem  lies  in  defining  a  physically  meaningful  elec¬ 
tron  distribution,  one  which  presumably  can  supply  a  constant  source  of 
free  energy  to  the  Kelvin-Helmholtz  instability,  and  which  is  plausible  in 
regard  to  the  complex  collisional  and  ionization  mechanisms  suspected  of 
generating  the  initial  magnetron  discharge.  This,  along  with  particle  sim¬ 
ulations  in  the  dynamic  regime,  will  be  the  subject  of  future  work. 
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Figure  3:  Schematic  of  the  planar  magnetron  discharge,  suggesting  gener 
at  ion  mechanisms. 


Figure  5:  Potential  profile  for  an  LS2  simulation  of  the  planar  magnetron, 
just  at  saturation,  K  =  Ke  =  4/9.  The  parameters  are  given  in  the  text. 
Full  line:  simulation  result,  Dots:  analytic  result,  Eqs.(15). 


Figure  6:  Ele< 
netron,  with  U 
Dots:  analytic 


Figure  7:  Current  versus  voltage  with  the  parameters  of  the  simulation  of 
Fig. (5).  We  consider  two  ion  temperatures,  with  nl0  adjusted  to  keep  the 
injection  rate  constant,  at  K  -  Ke  =  4/9. 


a 

gure  8:  Current  versus  magnetic  field,  for  the  same  parameters  as  in  fj 

g.(5).  The  value  of  the  critical  field,  Bc  =  4.0,  is  found  from  Eq.(l7).  $ 
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Figure  9:  Ion  scatter  plots,  vx,  versus  it,  for  varying  magnetic  fields.  The 
critical  field  is  Bc  =  4.0,  the  other  parameters  are  the  same  as  in  Fig. (5); 
(a)  B  =  0,  (b)  B  =  3.5,  (c)  B  =  4.5. 


SECTION  m:  CODE  DEVELOPMENT  AND  SOFTWARE  DISTRIBUTION 


Codes  which  we  haw  developed  a  modified  here  are,  generally  speaking,  available  for  the  ask¬ 
ing.  We  are  making  arrangements  with  our  Industrial  Liaison  Program  to  handle  distribution,  at 
ccst,  ctf  listing?,  tapes  arad/ar  diskettes  In  a  later  report  we  will  provide  a  list  erf  codes  and  related 
material  available  and  bow  to  obtain  such.  If  you  use  those,  we  simply  ask  that  yew  acknowledge  the 
aource  and  also  send  us  your  results  as  publi&ed  in  reports  or  journals. 

(Not  all  ctf  the  codes  mentioned  in  our  reports  are  available  for  distribution;  same  are  angle- 
purpose,  angle-use,  too  stylized  far  outside  use.) 


ES2:  An  Electrostatic  2d  Bounded  Simulation  Code  with  Injection 
and  External  Circuit  -  K.  Theilhaber. 

I.  INTRODUCTION. 

The  particle  simulation  code  ES2  has  been  designed  as  a  two-dimensional 
generalization  of  the  one-dimensional,  bounded  electrostatic  code  PDWl(l), 
and  imitates  or  is  inspired  by  many  features  of  this  one-dimensional  code.  It 
consists  of  a  bounded  simulation  region  (see  Fig .  ( 1 ) ) ,  containing  the  plasma 
particles  (electrons  and  ions),  with  periodic  boundary  conditions  on  the  parti¬ 
cles  and  fields  in  one  direction  (y),  and  finite  bounding  surfaces  in  the  other, 
perpendicular  direction  (x),  At  the  bounding  surfaces,  which  can  represent 
walls  or  plates  in  a  physical  device,  or  merely  open  boundaries  in  the  simula¬ 
tion  of  an  infinite  plasma,  the  code  allows  for  the  emission  and  absorption  of 
particles.  The  code  can  also  model  the  presence  of  an  external  circuit,  so  as  to 
permit  the  external  flow  of  charge  from  one  boundary  to  the  other,  subject  to 
the  desired  circuit  characteristics.  The  plasma  simulation,  while  incorporating 
only  electrostatic  self-forces,  also  allows  for  an  external  imposed  magnetic  field 
(Bo  in  the  figure),  which  can  be  set  at  an  arbitrary  angle  in  the  xyz  plane. 

ES2  has  been  specifically  designed  for  extending  the  areas  of  research  un¬ 
dertaken  with  PDWl,  in  particular  the  study  of  the  physics  and  dynamics  of 
the  sheathes  and  non-neutral  structures  which  arise  in  plasma-wall  interactions. 
In  the  MFE  environment,  it  can  run  on  both  Cray-l’s  (C  and  D  Machines),  on 
the  Cray  X-MP  (E  Machine),  and  is  being  adapted  to  run  on  the  Cray-2  (B 
Machine).  In  what  follows,  we  shall  first  give  a  brief  synopsis  of  the  normaliza¬ 
tions  introduced  in  the  physical  model.  We  shall  then  describe  the  numerical 
scheme  used  in  the  code,  and  the  various  options  available.  Finally,  we  give  an 
illustrative  example,  and  note  future  directions  of  development. 

n.  PHYSICAL  NORMALIZATIONS. 

In  the  electrostatic  approximation  assumed  by  ES2,  the  plasma  is  entirely 
described  by  the  coupled  Lorentz  and  Poisson  equations: 

=  —  (E  +  v  x  B)  ,  (1) 

at  m# 

V2<t>  =  -  — (n,  -  n,), 

<o 

E  =  -  V$s, 


(2) 

(3) 


where  a  is  the  species  index,  a  =  e,t,  and  where  n(  and  n,  denote  the  particle 
densities  (number  per  unit  volume).  We  assume  singly-charged  ions,  so  that 
9.  =  -<7«  =  «• 

In  our  simulation  model,  we  choose  e/me  =  1,  to  =  1  and  we  define  n, 
and  ne  as  the  number  of  particles  per  unit  area.  Thus,  the  equations  of  motion 
become: 


-(E  +  v  XB), 

(4) 

d  frig 

_v  =  —  (E  +v  x  B) , 
dt  m, 

(5) 
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For  instance,  let  us  consider  a  homogeneous  simulation  plasma  where  electrons 
and  ions  have  an  equilibrium  density  of  n  (per  unit  area),  and  where  the  elec¬ 
trons  are  initialized  with  thermal  velocity  o(e,  the  ions  with  thermal  velocity 
utl,  with  ss  (m,/m,)1/!iii,.  Then  the  plasma  and  cyclotron  frequencies  are 
simply: 


=  3 ,  (10) 

Wc,  =  —  B,  (11) 

m. 


If  we  then  choose  vte  =  1  and  assume  n  is  given,  then  by  setting  e  =  l/n 
we  insure  that  wpe  =  1  and  hence  that  Att,  =  1  as  well.  In  other  words,  n 
automatically  becomes  the  number  of  particles  per  Debye  square. 


m.  NUMERICAL  TECHNIQUES. 


We  shall  briefly  describe  the  numerical  schemes  use  in  ES2. 

The  field  solver  is  of  a  mixed  type,  similar  to  the  one  described  by  Birdsall 
and  Langdon  in  (3):  Poisson’s  equation  is  solved  by  using  a  Fourier  represen¬ 
tation  of  the  equation  in  the  y  direction,  that  is  by  transforming  from  y  to 
ky  using  an  FFT  (Fast  Fourier  Transform),  and  by  solving  a  finite-difference 
scheme  in  the  x  direction  for  each  of  the  Fourier  components.  A  final  inversion, 
from  ky  to  y  space  using  the  FFT,  completes  the  solution  process.  The  finite- 
difference  equation  for  the  ’d.c.’  (ky  =  0)  Fourier  component  requires  a  special 
treatment,  as  it  must  account  for  the  surface  charge  that  accumulates  on  the 
boundary  surfaces.  This  treatment  is  in  fact  nearly  identical  to  the  formulation 
of  Poisson’s  equation  in  PDWlf1).  For  the  higher  Fourier  components,  the 
boundary  conditions  are  sero,  and  the  solution  of  the  finite  difference  equations 
is  slightly  simpler(3). 

The  numerical  particles  are  ’square”  in  shape,  and  use  a  straightforward 
area-weighting  for  the  gather  /scatter  operations,  as  described  in  Ref.(3).  When 
the  particles  hit  a  boundary  they  are  absorbed,  and  their  charge  is  ’instan¬ 
taneously”  dispersed  as  a  uniform  surface  charge  on  the  conductor’s  surface. 
Similarly,  when  a  particle  is  injected  from  a  surface,  it  leaves  a  ’hole’  of  op¬ 
posite  charge  behind,  and  the  hole’s  charge  contribution  is  equally  accounted 
for.  A  repacking  algorithm,  copied  from  PDWl,  periodically  discards  the  coor¬ 
dinates  and  velocities  of  particles  which  have  been  absorbed  by  a  surface,  and 
thus  keeps  the  total  number  of  particles  in  the  system  at  any  one  time  within 
reasonable  bounds  of  sise. 

The  external  circuit  uses  an  algorithm  similar  to  the  one  described  in  (1). 
It  can  model  a  voltage  source,  a  resistor  and  an  inductor  in  series,  or  more 
simply  a  short  or  an  open  circuit  across  the  plasma  region. 


IV.  OPTIONS. 

Various  options  exist  for  setting  up  the  physical  conditions  of  a  simulation. 
The  magnetic  field  can  be  chosen  at  an  arbitrary  angle  to  the  xyz  axes.  While  in 
the  present  version  of  ES2  the  magnetic  field  is  constant,  the  algorithm  which 
moves  the  particles  can  be  modified  for  a  non-uniform  magnetic  field.  The 
program  can  be  initialized  with  the  simulation  region  full  or  empty  of  particles, 
with  prescribed  stationary  or  drifting  Maxwellians.  One  can  separately  specify 
the  injection  rates  of  each  species  at  each  bounding  surface.  The  distribution 
functions  of  the  emitted  particles  are  Maxwellians,  with  specified  temperatures 
and  drifts.  One  available  option  assumes  that  the  simulation  region  bathes 
in  a  "plasma  reservoir”,  of  given  density,  temperature  and  drift  velocities;  the 
injection  rate  is  then  automatically  adjusted  in  accordance  with  the  reservoir 
parameters. 
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At  the  present  time,  ES2  incorporates  diagnostics  in  the  form  of  scat¬ 
ter  plots,  graphs  of  the  distribution  functions,  field  plots  of  the  potential  and 
electric  fields  (cross-sections  and  contours),  plots  of  Fourier  components  of  the 
fields,  and  time  histories  of  energies,  circuit  quantities,  etc.  A  more  specialised 
diagnostic  provides  time-averaged  snapshots  of  grid  quantities. 


V.  AN  EXAMPLE  OF  AN  ES2  RUN. 

As  an  example  of  a  run  of  the  ES2  code,  we  shall  consider  the  two  dimen¬ 
sional  simulation  of  instability  induced  by  a  drifting  population  of  Maxwellian 
electrons,  in  the  presence  of  a  cold  (X)  «C  Tt )  population  of  ions.  This  is  the 
regime  for  the  so-called  ion  acoustic  instability,  predicted  by  linear  plasma  the¬ 
ory.  Through  particle  simulations,  we  can  hope  to  gain  insight  in  the  nonlinear 
dynamics  of  the  process  (,4,8). 

The  system  dimensions,  illustrated  in  Fig.(l),  are  Lz  =  32,  Lv  =  32,  with  a 
grid  spacing  Ax  =  1.0  and  Ay  =  1.0.  Thus,  the  numerical  mesh  has  Nxx  Nv  = 
33  x  32  =  1056  points.  The  system  is  initially  filled  with  equal  populations  of 
electrons  and  ions  (about  5000  of  each),  with  densities  ne  ~  rij  =  5  per  unit  area, 
electron  charge  e  =  0.2  and  mass  ratio  mi/me  =  100.  The  thermal  velocities 
are  u,e  =  1.0  and  =  0.01,  corresponding  to  "temperatures*  Te  =  vt2cme  =  0.2 
and  Ti  =  0.02,  so  that  T,  C  Te.  Furthermore,  the  electrons  are  initially  set 
drifting  with  a  velocity  Vde  =  1.0.  With  these  normalisations,  we  find  that 
Wpe  =  1,  cjp,  =  0.1,  c,  =  (T,/m,)1/3  =  0.1,  Xde  -  1-  In  particular,  we  have 
Lx  =  Lv  —  32Ade.  Finally,  in  this  run  we  chose  the  plasma  to  be  unmagnetised, 
with  19=0. 

The  injection  rates  are  set  to  those  that  would  obtain  if  the  simulation 
region  bathed  in  a  reservoir  of  plasma,  with  exactly  the  same  parameters  as 
those  outlined  above.  Thus,  in  the  absence  of  collisional  and  collective  effects, 
the  system  should  be  in  a  steady  state,  with  injection  at  the  boundaries  exactly 
replenishing  the  particles  lost  by  absorption  at  these  same  boundaries.  In  fact, 
this  is  not  the  case:  first,  because  even  when  =  0,  there  is  a  slow  but 
measurable  heating  of  the  ions  due  to  numerical  collisonality,  and  second,  when 
v<ie  /  0,  turbulence  develops  in  the  plasma. 

Finally,  the  external  circuit  is  set  to  a  short-circuit  condition:  we  require 
that  d>(0,  y)  =  4>(L ,,  y)  =  0  at  all  times. 

In  Fig. (2)  we  display  the  ion  phase  space  (  vzvs  x  )  at  (  =  140  (ui,nt  «  14). 
This  scatter  plot  shows  all  ions  in  the  system.  A  large  ion  vortex  has  formed 
about  i  =,  with  a  velocity  trapping  width  of  order  Avtr  a  0.08  3>  vt,  =  0.01. 
The  ions  in  the  vortex  are  trapped  in  a  large-amplitude  wave,  which  propagates 
in  the  direction  of  electron  drift  with  a  velocity  comparable  to  the  ion  acoustic 
velocity  c,  =  0.1.  The  two  dimensional  potential  structure  of  this  wave  is  shown 
in  Fig. (3):  this  is  a  perspective  view  of  the  potential  4>(x,y),  averaged  over  a 


time  interval  T  =  1.5u“‘  so  as  to  suppress  the  high-frequency  (w  es  uipt  =  1) 
fluctuations  of  the  fields.  Fig. (3)  displays  three  features  which  are  expected  in 
this  situation  (,45).  The  first  is  a  broad  potential  dip,  which  traps  the  ions 
observed  Fig. (2).  The  second  is  a  positive  potential  step,  downstream  from  the 
ion  vortex,  which  is  formed  by  the  asymmetric  reflection  of  electrons  by  the 
potential  dip.  The  third  feature  is  a  low-frequency  (w  <  Wpi)  turbulence,  which 
eventually  dominates  the  ’coherent’  features  of  the  potential  profile  that  we 
have  just  just  discussed. 

VI.  CONCLUSION. 

Further  development  of  ES2,  especially  in  connection  with  its  diagnostics, 
will  be  strongly  stimulated  by  its  continued  use  in  a  number  of  research  projects. 
We  plan  to  issue  a  first  ’User’s  Manual’  in  the  near  future. 
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time  interval  T  =  l.Sw"1  10  u  to  auppreu  the  high-frequency  (w  m  upe  =  1) 
fluctuation*  of  the  field*.  Fig. (3)  displays  three  feature*  which  are  expected  in 
thia  aituation  (,45).  The  fir*t  is  a  broad  potential  dip,  which  traps  the  ion* 
observed  Fig. (2).  The  second  is  a  positive  potential  step,  downstream  from  the 
ion  vortex,  which  is  formed  by  the  asymmetric  reflection  of  electrons  by  the 
potential  dip.  The  third  feature  is  a  low-frequency  (w  <  w,*)  turbulence,  which 
eventually  dominates  the  "coherent*  features  of  the  potential  profile  that  we 
have  just  just  discussed. 

VI.  CONCLUSION. 

Further  development  of  ES2,  especially  in  connection  with  its  diagnostics, 
will  be  strongly  stimulated  by  its  continued  use  in  a  number  of  research  projects. 
We  plan  to  issue  a  first  ’User’s  Manual*  in  the  near  future. 
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B .  Energy  Conserving  Reflection  Algorithm 
Lou  Ann  Schwager 

I.  Introduction 

The  simulation  of  plasma  reflecting  off  a  boundary  requires  the  correct  formulation  of 
an  elastic  reflection  of  a  plasma  ion  off  a  wall  or  collector  plate.  Conservation  of  total  energy 
at  the  event  of  reflection  is  crucial  when,  in  a  steady  state  system,  particles  injected  at  a 
source,  traverse  a  potential  path,  perfectly  reflect  at  a  wall,  and  return  to  the  origination 
point.  In  this  model  these  particles  should  have  no  gain  or  loss  in  energy.  In  a  real  model, 
some  of  the  energy  of  the  reflected  particle  is  imparted  to  the  wall  in  inelastic  reflection. 
Hence  we  provide  this  exercise  as  a  reference  for  these  cases. 

During  one  time  step  of  the  simulation,  a  particle  incident  on  a  wall  passes  beyond 
it.  The  penetration  depth  depends  on  the  previous  location  and  present  velocity  of  the 
particle.  This  velocity  is  determined  with  a  leap-frog  mover  at  one-half  time  step  earlier. 
On  reflection  the  position  of  the  particle  is  mirrored  back  through  the  wall.  As  a  first 
guess,  we  give  the  new  velocity  the  same  magnitude  but  opposite  sign  of  the  old  velocity. 
This  assumption  does  not  conserve  energy  because  the  particle  has  been  reflected  across  a 
potential  difference.  The  following  derivation  provides  a  correction  in  this  simply  reflected 
velocity  which  includes  the  effect  of  the  potential  change  and  one-half  time  step  correction. 

II.  Model 

A  particle  of  mass  m,  charge  q,  and  velocity  u  moving  along  x  in  a  time-independent 
electric  field  E  has  constant  energy  £ , 


The  acceleration  a  on  the  particle  is 

-*(=)•  m 
We  assume  this  acceleration  is  constant  within  one  spatial  grid  Ax  of  the  wall  and  over  one 
time  step  At  prior  to  reflection. 

If  the  subscript  p  denotes  the  time  step  prior  to  wall  penetration,  then  the  simulation 
gives  us  xp  and  vp_  ^ .  At  the  next  time  step  when  the  particle  is  moved,  the  new  velocity 
and  position  become 

Vt-J  =  vp~i  +aAt  (3) 

and 

Xp+i  =  Xp  +  Vp+^A  t.  (4) 

If  Xp+ 1  exceeds  the  system  length  L  and  the  particle  is  reflected  (as  opposed  to  added  to 
the  wall  surface  charge),  then  the  new  position  and  velocity  are  xp+1  and  .  We  choose 
to  reposition  the  particle  from  the  wall  at  the  distance  it  would  have  penetrated  the  wall  1 
Hence, 

Zp+i  =  2  L  xp+ 1-  (5) 

This  scheme  is  shown  in  Fig.  1. 

We  next  correct  the  simple  velocity.  The  particle,  first  accelerated  by  E  at  the  wall 
to  Xp+i,  is  now  decelerated  by  the  same  E  to  £p+  j.  We  correct  the  simple  velocity  with  6. 
where 

+6.  (6) 

Equating  £p  with  £p+ i  determines  the  value  for  .  This  energy  balance  gives 


Second-order  approximations  for  v3  at  the  full  time  step  are 


tV  =  lWitV-i 


Vi=V!Vj 


The  next  velocity  depends  on  the  reflected  velocity  by 


ty+j  =  tyj  +  oAt. 


Finally  we  substitute  Eqs.  (3)-(6)  and  (8)-(10)  into  the  energy  balance  of  Eq.  (7).  The 


solution  for  6  from  this  energy  balance  is  in  the  quadratic  equation: 


0  =  62  +  6(aAt  -  2iy.j )  +  4o(ryi  -  L)  -  2aiy+j  At. 


To  evaluate  6 ,  we  consider  the  ordering  of  terms  in  Eq.  (11)  Because  the  leap-frog 
mover  estimates  position  and  velocity  through  a  first-order  central  difference,  these  ap¬ 
proximations  (Eqs.  (3),  (4),  and  (10))  are  of  0((At)3).2  Simple  reflected  velocity  is  0(1) 
because  as  At  —  0,  then  ry  s  — *  L  and  ty  j  — »  -ty  ^  Finally  the  v3  approximation  in 
Eqs.  (8)  and  (9)  introduces  an  error  with  0((At)2) ,  which  then  is  the  lowest-order  error  of 
Eq  (11).  This  determines  that  6  is  O(At).  If  we  neglect  terms  with  0((At)2)  or  greater. 


then  the  second-order  approximation  for  6  becomes 


6  =  - (x_+i  -  L)  -  aAt 


Note  that  for  no  wall  charge,  a  =  0,  and  so  ty  j  =  — ty  j 


III.  Conclusion 


The  maximum  wall  penetration,  Xp-n  —L  =  ty  ^  At,  fixes  6  at  +aAt  Minimum 
penetration,  Xp+i  =  L,  generates  the  smallest  correction  with  6=  —aAt.  Because  the  value 


WWW 


vw 


PV 


of  Xp+i  is  equally  probable  across  the  range  from  L  to  L+iy.^  At,  the  average  value  of  6 
after  many  reflections  is  zero.  On  the  other  hand,  the  correction  6  generates  a  spread  in 
velocities  of  the  reflected  ions.  Consequently,  if  a  substantial  electric  field  develops  at  the 
wall  and  we  use  only  simple  reflection,  the  reflected  ion  flux  would  be  erroneously  cooled. 

Typically  for  a  low  temperature  plasma,  the  wall  absorbs  about  20%  of  the  incident 
particle  energy.  The  specific  fraction  depends  on  the  mass  ratio  of  the  incident  particle 
and  wall  target.3  For  this  example,  the  velocity  of  the  reflected  particle  is  about  90%  of  its 
incident  velocity.  Specifically  in  the  simulation  the  new  reflected  velocity  would  be 

=  -°-9  V4  +  6-  (13) 
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Direct  Implicit  Particle  Simulation  of  Tandem 
Mirrors.*  B.  I.  COHEN,  M.  E.  STEWART,  and  R.  P. 
FREIS,  Lawrence  Livermore  National  Laboratory,  and  L.  A. 
STRUGALA  and  C.  K.  BIRDSALL,  University  of  California, 
Berkeley  -  TESS  is  a  one-dimensional  axial  particle  code  for 
simulating  thermal-barrier  formation  and  loss  of  plugging  in 
a  tandem  mirror  plasma.  TESS  calculates  the  self-consistent 
ambi polar  potential  from  data  collected  from  advancing  par¬ 
ticle  ion  and  electron  trajectories  from  the  midplane  of  the 
center  cell  to  the  end  of  the  device.  By  performing  a  direct  im¬ 
plicit  solution  for  the  electric  potential,  we  are  allowed  to  use 
a  long  time  step.  An  accurate  representation  of  passing  and 
trapped,  relativistic  electron  and  ion  guiding-center  dynam¬ 
ics  is  rendered  by  conventional  particle  simulation  and  Monte 
Carlo  techniques.  Charge  exchange  and  ionization  processes 
are  represented  probabilistically:  The  wave  heating  effects  of 
ECRH  and  ICRF  are  modeled  with  the  quasilinear  formalism 
of  Bernstein  and  Baxter.  Transit  loss  and  self-consistent  bi¬ 
nary  Coulomb  collisions  that  conserve  energy  and  momentum 
are  included.  A  probability  for  radial  loss  is  calculated  analyt¬ 
ically  for  classical  and  neoclassical  radial  transport. 

Work  performed  for  U.S.D.O.E.  by  Lawrence  Livermore  National 

Laboratory  under  contract  W-7405-ENG-48. 


Simulations  of  the  Beam-Plasma  Instability 
In  the  Turbulent-Trapping  Regime,  K.  THE1LHABER,  U.C. 
Berkeley.  G.  LAVAL  and  D.  PESME,  Ecole  Poly  technique. 
France .  —  In  the  weak  warm-beam-plasma  instabilitv, 
the  "turbulent- trapping"  regime ^  is  defined  byys(jJl^ 
>  where  &  is  the  growth-rate  and  U  the  resonance- 
broadening  frequency.  In  this  regime,  trapping  domi¬ 
nates  the  phase-space  structure  of  the  distribution 
function,  leading  to  the  generation  of  a  harmonic  se¬ 
ries  of  quasimodes  resonant  with  the  beam  particles, 
and  to  a  theoretical  increase  of  growth-rates  and 
diffusion  above  the  quasilinear  values.  While  our 
simulations  confirm  the  detailed  structure  of  the 
correlation  functions  predicted  by  the  turbulent- 
trapping  theory,  they  do  not  show  the  expected  enhan¬ 
cement  of  the  growth-rates.  We  shall  discuss  these 
results,  and  link  them  to  a  numerical  study  of  velo¬ 
city-space  diffusion  in  the  turbulent  fields.  This 
will  lead  to  a  revision  of  the  Fokker- Planck  equation 
for  the  correlation  function. 

1.  G.  Laval  and  D.  Pesme,  Pbys.  Rev.  Lett.,  53,3,  270 
(1984).  — 


Supported  by  DOE  Contract  DE-AT03-76ET5 >0b4  and  ONR 
Contract  N000 1 4-77-C-0578  . 


Single  Ended  Q  Machine  Equilibria  -  Experiment 
vs.  Theory.*  M.  OERTL,  P.  KRUMM,  and  S.  KUHN,  Inns¬ 
bruck  University,  Austria,' — Equilibrium  studies  of  the 
single  ended  Q  machine  with  negatively  biased  cold 
plate,  based  on  a  collisionless  plane  diode  model, 
have  predicted1  that,  for  increasing  q  (neutral  flux 
towards  the  hot  plate),  both  ne  (plasma  density)  and 
(ion  current  density)  should  first  increase  mono- 
tonically  (regime  1',  monotonically  decreasing  poten¬ 
tial  profile);  after  reaching  some  "transition  point", 
ne  should  decrease  slightly  while  should  remain  al¬ 
most  constant  (regime  2',  potential  profile  with  maxi¬ 
mum  near  the  hot  plate).  In  the  present  experiment, 
this  prediction  has  been  checked  by  monitoring  ne  and 
j £  while  increasing  q,  for  differnt  values  of  the  hot 
plate  temperafure  T.  While  regime  1'  and  the  transit¬ 
ion  point  are  recovered  in  accordance  with  theory,  re¬ 
gime  2'  shows  strongly  increasing  and  slightly  in¬ 
creasing  ne.  This  evidence,  along  with  recent  simula¬ 
tion  results, ^  suggests  a  slight  modification  of  the 
theory  so  as  to  allow  for  trapped  electrons  in  the  re¬ 
gion  of  the  potential  maximum  near  the  hot  plate. 

*Work  supported  by  Austrian  Research  Funds  Contracts 

S-18/02 ,  S-18/03,  and  PS178. 

is.  Kuhn,  Plasma  Phys.  881  (1981) 

2S..  Kuhn  et  al..  Proceedings  of  the  1984  Int.  Conf.  on 
Plasma  Physics,  Lausanne,  Contributed  Papers,  p.  104. 


Equilibrium  Corrections  for  Electron  Hole 
Filling:  Simulation  vs.  Theory  ,  P .  C .  GRAY ,  T .  L . 
CRYSTAL,  S.  KUHN*  and  C.K.  BIRDS ALL,  U.C.  Berkeley-- 
A  recent  quantitative  theory  for  axial  one  dimensional 
equilibria  of  a  collisionless  single-ended  Q-machine 
predicts  the  maximum  and  plasma  potential  for  the  case 
when  the  considered  particle  source  is  taken  ion  rich 
and  the  over  all  bias  is  held  negative1.  In  this  case 
a  single  positive  potential  extremum  is  near  the 
cathode  and  the  steady  state  theory  assumes  an 
electron  hole  in  phase  space  from  which  all  electrons 
are  forbidden.  However,  particle  simulations  of  this 
bounded-plasma  configuration  recover  very  different 
values  of  potential,  and  the  assumed  electron  hole  is 
found  to  be  filled  in2:  the  electron  distribution 
f(v,x)  is  measured  to  be  Maxwell-Boltzmann,  f($,x). 
This  result  may  also  explain  some  recent  equilibrium 
measurements2. —  Work  supported  by  DOE  Contract 
DE-AS03-76-FOOO 34-DE-AT03- 76ET53064 ,  ONR  Contract 
N000 1 4-77-C-O  578  and  Austrian  Research  Funds. 
Simulations  performed  at  NMFECC. 

1.  Kuhr,  S. .Plasma  Phys.,  Vol.  23,  pp.  881-902,  1981. 

2.  Kuhn,  S.,  Prbc.  I.C.P.P.,  Lausanne  Swi . ,  Vol.  1, 
pp .  104,.  'June  1984. 

3.  Ocrtl,  M.  et  al . ,  this  conference. 
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A  Slroulation  Study  of  the  Dependence  of  the 
Plasma  Sheath  Impedance  on  Frequency ,  C.  K.  BIRDSALL, 

P.  C.  CRAY  and  T.  L.  CRYSTAL,  U.  C.  Berkeley-- 
A  one  dimensional  electrostatic  particle  code  bounded 
with  an  external  circuit^has  been  used  to  model  the 
parallel  sheath.  Static  potential  profiles  agree  with 
simple  Boltzmann  electron  theory.  The  model  has  been 
used  to  measure  the  frequency  dependence  of  the  sheath- 
impedance  about  the  floating  potential,  for  frequencies 
0<  u  <  5uipe  .  The  impedance  includes  effects  from 
potential  drops  at  both  the  plasma  source  and  the 
collector.  For  a  short  system  at  low  frequency  the 
impedance  appears  resistive,  with  Rjq  =  1/g^Q  where 
g^Q  is  obtained  from  dc  measurement  of  t)J/d^  at  the 
floating  potential.  A  parallel  resonance  is  found  at 
u>“  u;pe  of  the  plasma  region,  due  to  the  vacuum  C  of 
bounded  system  and  L  associated  with  electron  inertia. 
It  is  necessary  to  shunt  the  sheath  resistance  with 
a  large  capacitance  in  order  to  account  for  the  ion 
imobility  at  w  £  wpi •  The  overall  result  is  an 
approximate  RLC  equivalent  circiut  for  the  sheath. 

For  a  longer  system,  $2(x)  for  various  choices  of 
driving  frequency  has  been  computed  showing  the  spatial 
behavior  of  the  signal. — Work  supported  in  part  by  DOE 
contract  DE-AS03-76-F00034-DE-AT03-76ET53064 ,  ONR 
contract  N00014-77-C-0578  and  a  gift  from  Varian 
Associates . 

1.  Wm.  Lawson,  U.C.  Berkeley,  Memorandum  No.  UCB/ERL 
M84/45,  June  13,  1984. 


A  Model  of  the  Plasma-Sheath  Region 
Including  Secondary  Electron  Emission  and  Ion  Re- 
f lection .  L.A.  Schwager,  C.K.  Birdsall,  Electronics 
Research  Lab,  U.C.  Berkeley.* — The  interaction  of  a 
medium  temperature  plasma  with  a  collector  plate  is 
studied  allowing  emission  of  secondary  electrons  and 
reflection  of  plasma  ions.  These  effects  alter  the 
electrostatic  potential  drop  across  the  sheath  re¬ 
gion  and  plasma  transport  to  the  plate.  The  compu¬ 
ter  code  PDW1  models  this  behavior,  solving  for 
particle  motion- in  the  self-consistent  electrostatic 
field,  in  one-dimension.  Particle  ions  and  elec¬ 
trons,  which  are  injected  continuously  from  a  hot 
plate  into  a  plasma  region,  flow  to  a  floating  plate 
which  reflects  ions  or  emits  secondary  electrons. 

Our  ion  reflection  simulation  results  compare  well 
with  our  static  cold  ion  theory.  We  also  observe 
an  ion-ion  two-streaming  instability  as  predicted 
by  our  time-dependent  theory.  Our  secondary  elec¬ 
tron  results  compare  well  with  previous  theory  and 
demonstrate  charge  saturation  as  predicted. ^ 
^Supported  by  DOE  contract  DE-AT03-76ET530A4 . 

^G.D.  Hobbs  and  J.A.  Wesson,  Plasma  Physics,  Vol.9, 

(  1967)  . 
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Experimental  and  Simulation  Study  of  Planar 
Magnetron  Discharge**  A.  WENDT,  P.C.  GRAY,  H.  MEUTH, 
M.A.  LIQJERMAN  and  C.K.  BIRDSALL,  U.  C.  Berkeley— 

An  apparatus  has  been  constructed  for  Investigation 
of  plasma  processes  In  planar  magnetron  discharges. 
Objectives  Include  finding  the  mechanisms  and  values 
of  ion  acceleration  and  electron  current,  aa  well  as 
the  role  of  waves  and  fluctuations.  The  experiment 
consists  of  independently  positioned,  moveable  1 

parallel  plate  electrodes  in  a  cylindrical  chamber 
and  diagnostics  that  allow  for  axially  as  well  as 
radially  resolved  measurements  of  temperature,  density 
and  fluctuations.  We  are  also  using  an  electrostatic 
bounded  particle  code  (1  space,  3  velocity  dimensions) 
to  simulate  the  magnetron  discharge.  The  code  >  1 

includes  modeling  of  bulk  ionization  and  electron- 
neutral  collisions.  Comparision  of  preliminary 
experimental  measurements  and  simulation  results  will 

be  given. 

* 

Research  supported  by  a  gift  from  Varian  Associates. 


A  Simulation  Study  of  Effects  of  Emission 
Noise  in  a  Thermionic  Diode  ,  Wm.  LAWSON,  P.  C.  GRAY, 
E.  A.  ADLER*and  C.  K .  BIRDSALL,  U.  C.  Berkeley  -- 
The  transition  between  the  space-charge-limi ted  and 
temperature-limited  operating  regimes  of  a  thermionic 
diode  is-  not  well  understood.  It  has  been  found 
experimentally  that  the  current-voltage-temperature 
characteristic  fits  a  simple  universal  formula*. 

The  wide  spread  applicability  of  this  formula  suggests 
that  the  physics  of  the  transition  involves  some 
universal  mechanism;  it  is  not  solely  a  question  of 
the  surface  properties  of  the  emitter.  It  has  been 
suggested  that  the  interaction  of  emission  noise  with 
the  space-charge  cloud  just  outside  the  cathode 
surface  may  affect  the  observed  dc  current-voltage 
characteristic  near  the  transition.  This  effect  has 
been  studied  using  the  PDW1  simulation  code^.  We 
will  discuss  the  effect  of  varying  levels  of  emission 
noise  on  both  the  current-voltage  characteristics  and 
the  dynamics  of  the  potential  minimum. — Work  supported 
in  part  by  Hughes  Aircraft  Company. 

1.  R.T.  Longo,  I EDM  1978,  p . 1 5 2 . 

2.  Wm.  Lawson,  II.  C.  Berkeley,  Memorandum  No.  UCB/ERL 
M84/45,  June  13,  1984. 

♦Hughes  Aircraft  Coapany,  Torrance. 
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